#ifndef __ARRAY3D_H__
  #define __ARRAY3D_H__

#include <string>
#include <sstream>

#include "exceptions.h"
#include "global_variables.h"

using namespace std;

//! Three dimensional array

/*!
A simple array implementation.
For arrays larger than 1GB some changes may be necessary.

Implementation date: 2008

Usage:

\code
array3d<int> a(2, 3, 4, 0); // with 2 pages, 3 rows and 4 columns all initialized with 0

a(0, 0, 0) = 17;

a(0, 1, 0) = 21;

a(1, 0, 3) = -4;
\endcode
*/

template <class T>
class array3d {
private:
  T *content;
  unsigned int pages;
  unsigned int rows;
  unsigned int cols;

  void allocate();
  void release();
  void hard_copy(const array3d &a);

public:
  //! resets the array to an empty array and releases allocated memory
  void clear();

  //! initializes all elements of the array to init_value
  void initialize(const T &init_value);

  //! returns the number of pages
  unsigned int get_page_count() const;

  //! returns the number of rows
  unsigned int get_row_count() const;

  //! returns the number of columns
  unsigned int get_column_count() const;

  //! reconstructs the array and deletes its content
  void reconstruct(unsigned int pages, unsigned int rows, unsigned int cols);

  //! resizes the array and keeps its content intact
  void resize(unsigned int pages, unsigned int rows, unsigned int cols, const T &init_value);

  //! index operator (may throw an exception)
  T &operator()(unsigned int pages, unsigned int row, unsigned int col) const;

  //! assignment operator
  array3d &operator=(const array3d &a);

  //! copy constructor
  array3d(const array3d &a);

  //! default constructor: an empty array, reconstruct or resize must be called to enlarge the array before it can be used
  array3d();

  //! constructor: does not initialize the array content
  explicit array3d(unsigned int pages, unsigned int rows, unsigned int cols);

  //! constructor: initializes the array content
  explicit array3d(unsigned int pages, unsigned int rows, unsigned int cols, const T &init_value);

  //! destructor (non-virtual: please do not inherit from this class)
  ~array3d();
}; // class array3d

//////////////////////////////////////////////////////////////////////////////

template <class T>
void array3d<T>::allocate() {
  release();
  content = new T[pages * rows * cols];
}

template <class T>
void array3d<T>::release() {
  delete[] content;
  content = 0;
}

template <class T>
void array3d<T>::hard_copy(const array3d &a) {
  this->pages = a.pages;
  this->rows = a.rows;
  this->cols = a.cols;

  allocate();
  const unsigned int total = pages * rows * cols;

  #pragma omp parallel for if (total > MINIMUM_SIZE_LIMIT_OMP)
  for (unsigned int i = 0; i < total; ++i) // do not use memcpy (it may cause problems with non-standard datatypes)
    this->content[i] = a.content[i];
}

template <class T>
void array3d<T>::clear() {
  release();
  rows = 0;
  cols = 0;
  pages = 0;
}

template <class T>
void array3d<T>::initialize(const T &init_value) {
  const unsigned int total = pages * rows * cols;

  #pragma omp parallel for if (total > MINIMUM_SIZE_LIMIT_OMP)
  for (unsigned int i = 0 ; i < total; ++i)
    content[i] = init_value;
}

template <class T>
unsigned int array3d<T>::get_page_count() const {
  return pages;
}

template <class T>
unsigned int array3d<T>::get_row_count() const {
  return rows;
}

template <class T>
unsigned int array3d<T>::get_column_count() const {
  return cols;
}

template <class T>
void array3d<T>::reconstruct(unsigned int pages, unsigned int rows, unsigned int cols) {
  this->pages = pages;
  this->rows = rows;
  this->cols = cols;

  allocate();
}

template <class T>
void array3d<T>::resize(unsigned int pages, unsigned int rows, unsigned int cols, const T &init_value) {
  if (this->pages >= pages && this->rows >= rows && this->cols >= cols)
    return; // do not shrink the array now!

  const unsigned int total = pages * rows * cols;
  T *tmp = new T[total];

  #pragma omp parallel for if ((pages * rows) > MINIMUM_SIZE_LIMIT_OMP) collapse(2)
  for (unsigned int i = 0 ; i < pages; ++i) {

    for (unsigned int j = 0 ; j < rows; ++j) {
      for (unsigned int k = 0 ; k < cols; ++k) {
        if (i < this->pages && j < this->rows && k < this->cols)
          tmp[i * (rows * cols) + (j * cols) + k] = content[i * (this->rows * this->cols) + j * this->cols + k];
        else
          tmp[i * (rows * cols) + (j * cols) + k] = init_value;
      } // for k
    } // for j

  } // for i

  this->pages = pages;
  this->rows = rows;
  this->cols = cols;
  release();
  content = tmp;
}

template <class T>
T &array3d<T>::operator()(unsigned int page, unsigned int row, unsigned int col) const {
  if (page >= this->pages || row >= this->rows || col >= this->cols)
    throw index_exception(global_empty_string);

  return content[page * (this->rows * this->cols) + (row * this->cols) + col];
}

template <class T>
array3d<T> &array3d<T>::operator=(const array3d &a) {
  if (&a == this)
    return *this;

  hard_copy(a); // content will be released first
  return *this;
}

template <class T>
array3d<T>::array3d(const array3d &a) : content(0), pages(a.pages), rows(a.rows), cols(a.cols) {
  hard_copy(a);
}

template <class T>
array3d<T>::array3d() : content(0), pages(0), rows(0), cols(0) {
  // nothing to do
}

template <class T>
array3d<T>::array3d(unsigned int pages, unsigned int rows, unsigned int cols) : content(0), pages(pages), rows(rows), cols(cols) {
  allocate();
}

template <class T>
array3d<T>::array3d(unsigned int pages, unsigned int rows, unsigned int cols, const T &init_value) : content(0), pages(pages), rows(rows), cols(cols) {
  allocate();
  initialize(init_value);
}

template <class T>
array3d<T>::~array3d() {
  release();
}

#endif
